Read in data analysis results
### Read in data analysis results
data_folder = paste0(params$data_folder, params$subj_name, "/")
path_vec = list.files(data_folder, full.names = TRUE, recursive = TRUE)
res_list = list()
for(m in 1:1){
path = path_vec[m]
### Read in results from data
load(path)
res_list[[m]] = res
}
### Collect from all subjects: clusters_list, center_pdf_array, v_vec
clusters_list = lapply(res_list, function(res) res$clusters_list)
center_pdf_array_list = lapply(res_list, function(res) res$center_pdf_array)
v_vec_list = lapply(res_list, function(res) res$v_vec)
t_vec = res$t_vec
n0_mat_list = lapply(v_vec_list, function(v_vec) n0_vec2mat(n0_vec = v_vec/(t_vec[2]-t_vec[1])))
# ### Fix spike intensity
# for (ind in 1:length(center_pdf_array_list)) {
# if (dim(center_pdf_array_list[[ind]])[3]<length(t_vec)) {
# center_pdf_array_list[[ind]] = get_center_pdf_array_v2(edge_time_mat_list = edge_time_mat_list[ind], clusters_list = clusters_list[ind], n0_mat_list = n0_mat_list[ind], t_vec = t_vec)
# }
# }
Visualize estimated intensities
if (!exists("t_vec")) {
t_vec = seq(0,336,1)
}
### Get connecting probabilities
order_clus_1 = order(rowSums(apply(center_pdf_array_list[[1]],
MARGIN = c(1,2),
function(vec)sum(vec*(t_vec[2]-t_vec[1])))),
decreasing = TRUE)
### Visualize estimated connecting intensities
clus_size_list = lapply(clusters_list, function(clusters)sapply(clusters,length))
g1 = plot_pdf_array_v2(pdf_array_list = center_pdf_array_list[[1]][order_clus_1,order_clus_1, ,drop=FALSE],
pdf_true_array = center_pdf_array_list[[1]][order_clus_1,order_clus_1, ,drop=FALSE],
clus_size_vec = clus_size_list[[1]][order_clus_1],
t_vec = t_vec, x_lim = c(), y_lim = c(0,max(unlist(center_pdf_array_list[[1]]))))
g1
Visualize estimated time shifts
### Get memberships for ALL nodes
membership_list = lapply(clusters_list, clus2mem)
membership_list[[1]] = clus2mem(clusters_list[[1]][order_clus_1])
# for (m in 1:length(membership_list)) {
# N_node = nrow(locs_mat_list[[m]])
# mem_tmp = rep(0,N_node)
# mem_tmp[non_iso_inds_list[[m]]] = membership_list[[m]]
# membership_list[[m]] = mem_tmp
# }
### Plot time shifts distribution
# data = tibble(v_vec=unlist(v_vec_list),
# membership=c(membership_list[[1]][non_iso_inds_list[[1]]], membership_list[[2]][non_iso_inds_list[[2]]]),
# subj=rep(1:length(v_vec_list), times=sapply(v_vec_list,length)))
data = tibble(v_vec=unlist(v_vec_list),
membership=c(membership_list[[1]]),
subj=rep(1:length(v_vec_list), times=sapply(v_vec_list,length)))
g4 <- data %>%
ggplot(aes(x=v_vec,group=as.factor(membership),
color=as.factor(membership)))+
geom_density( aes(), alpha=0.6, position = 'dodge') +
scale_fill_manual(values=palette()[1+1:3]) +
theme_bw()+
theme(legend.position = "none") +
# facet_wrap(~subj) +
scale_x_continuous() +
theme(axis.ticks.y = element_blank(), axis.text.y = element_blank(), axis.title.y = element_blank()) +
xlab('Time shifts (min)')
print(g4)
## Warning: Width not defined. Set with `position_dodge(width = ?)`
### Get Wan's result
dat=readMat(paste('../../Raw_data/to_shizhe/','Leader_',substr(params$subj_name, start=nchar(params$subj_name)-7, stop=nchar(params$subj_name)),'.mat',sep=''))
activeTime_vec = dat$activeTime * 60
patternTime_vec = dat$patternTime * 60
# ### Let time shift for non-active nodes be NA
# for (m in 1:length(v_vec_list)) {
# N_node = nrow(locs_mat_list[[m]])
# v_vec_tmp = rep(NA,N_node)
# v_vec_tmp[non_iso_inds_list[[m]]] = v_vec_list[[m]]
# v_vec_list[[m]] = v_vec_tmp
# }
ggplot(data = data.frame(activeTime_vec=activeTime_vec,
v_vec=unlist(v_vec_list),
mem=as_factor(unlist(membership_list))
)) +
geom_point(aes(x=activeTime_vec, y=v_vec, color=mem)) +
# scale_color_manual(values=palette()[c(8,2:6)]) +
scale_color_manual(values=palette()[c(2:6)]) +
xlab('Activation time by Wan (pp. e7)') +
ylab('Our estimation')
## Warning: Removed 10 rows containing missing values (geom_point).
ggplot(data = data.frame(activeTime_vec=activeTime_vec,
patternTime_vec=patternTime_vec,
v_vec=unlist(v_vec_list),
mem=as_factor(unlist(membership_list))
)) +
geom_point(aes(x=patternTime_vec, y=v_vec, color=mem)) +
# scale_color_manual(values=palette()[c(8,2:6)]) +
scale_color_manual(values=palette()[c(2:6)]) +
xlab('Patterning time by Wan (pp. e7)') +
ylab('Our estimation')
## Warning: Removed 18 rows containing missing values (geom_point).
{r} # ind_vec = which(unlist(v_vec_list)!=0 & activeTime_vec==1) # ind = ind_vec[1] # # for (i in 1:min(length(ind_vec),3)) { # ind = ind_vec[i] # data_trace = tibble(dFF=dFF_mat[ind,],time=(1:ncol(dFF_mat)*0.25)/60) # g = ggplot(data_trace, aes(x=time, y=dFF)) + # geom_line(alpha=0.3) + # # xlim(unlist(v_vec_list)[ind]+c(-1,1)) + # geom_vline(xintercept = unlist(v_vec_list)[ind]) # print(g) # } # # #
Visualize estimated clusters with spatial location
### Spatial location ----
data = as.data.frame(reduce(locs_mat_list,rbind))
colnames(data) = c('AP','LR','DV')
data = cbind(data,
membership=as.factor(unlist(membership_list)),
subj=rep(1:length(membership_list),
times=sapply(membership_list,length)))
# lims = data.frame(coord=c("AP","DV","LR"), xmin=c(-0.5, 0, -2),xmax=c(9,140,3))
g<- data %>%
pivot_longer(cols = c('AP','LR','DV'), names_to = "coord", values_to = "loc") %>%
ggplot(aes(x=loc,
group=membership,
color=membership))+
geom_density( aes(), alpha=0.6, position = 'dodge') +
scale_color_manual(values=palette()[c(8,2:6)]) +
theme_bw()+
theme(legend.position = "none") +
facet_wrap(~coord, scales='free',ncol=1) +
xlab('') +
# scale_x_continuous(limits = c(-2, 6))+
theme(axis.ticks.y = element_blank(), axis.text.y = element_blank(), axis.title.y = element_blank()) +
scale_y_continuous(position = "right")
print(g)
## Warning: Width not defined. Set with `position_dodge(width = ?)`
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:igraph':
##
## groups
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
fig = plot_ly(data, x = ~LR, y = ~AP, z = ~DV, color = ~membership,
# colors = c(palette()[c(8,2:6)]),
colors = c(`0`=palette()[8],
`1`=palette()[2],
`2`=palette()[3],
`3`=palette()[4],
`4`=palette()[5],
`5`=palette()[6]),
type = "scatter3d")
fig
## No scatter3d mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
Visualize estimated clusters with mnx+/- indicator
### mnx decomposition ----
data = tibble(subj=rep(1:length(membership_list),sapply(membership_list,length)),
membership=unlist(membership_list),
mnx=as.factor(unlist(mnx_vec_list)),
islet=NA)
data_2 = data %>%
mutate(mnx=factor(mnx,levels=c(1,0)), islet=as.factor(as.numeric(islet)),
membership=factor(membership, levels=c(0,max(membership):1)))
data_2_left = data_2 %>% filter(subj==1) %>% mutate(membership=fct_drop(membership))
data_2_right = data_2 %>% filter(subj==2) %>% mutate(membership=fct_drop(membership))
fill_color = c(`0`=palette()[8],
`1`=palette()[2],
`2`=palette()[3],
`3`=palette()[4],
`4`=palette()[5],
`5`=palette()[6])
### Left spine
g1 = data_2_left %>%
ggplot(mapping = aes(fill=membership,x=mnx)) +
scale_fill_manual(values=fill_color) +
geom_bar(position="fill") +
geom_text(data = data_2_left %>%
count(mnx, membership) %>% group_by(mnx) %>%
summarise(perc=round(n/sum(n),3), membership=membership) %>%
arrange(desc(membership), .by_group = TRUE) ,
aes(x = mnx, label = scales::percent(perc), y=perc, fill=NULL),
position = position_fill(0.5), size=2.5)+
scale_x_discrete(labels=c("mnx+","mnx-"))+
xlab("")+
theme_bw()+
theme(legend.position = 'none',
panel.grid = element_blank(),
# axis.text = element_blank(),
axis.text.y = element_blank(),
axis.title = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank())
## `summarise()` has grouped output by 'mnx'. You can override using the `.groups` argument.
g1
